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Abstract 

Seeing the Earth crust as crisscrossed by fauhs filled with fluid at close to lithostatic pres- 
sures, we develop a model in which its elastic modulii are different in net tension versus com- 
pression. In constrast with standard nonlinear effects, this "threshold nonlinearity" is non- 
perturbative and occurs for infinitesimal perturbations around the lithostatic pressure taken as 
the reference. For a given earthquake source, such nonlinear elasticity is shown to (i) rotate, 
widen or narrow the different lobes of stress transfer, (ii) to modify the 2D-decay of elastic 
stress Green functions into the generalized power law where 7 depends on the azimuth 
and on the amplitude of the modulii asymmetry. Using reasonable estimates, this implies an en- 
hancement of the range of interaction between earthquakes by a factor up to 5 — 10 at distances 
of several tens of rupture length. This may explain certain long-range earthquake triggering and 
hydrological anomalies in wells and suggest to revisit the standard stress transfer calculations 
which use linear elasticity. We also show that the standard double-couple of forces representing 
an earthquake source leads to an opening of the corresponding fault plane, which suggests a 
mechanism for the non-zero isotropic component of the seismic moment tensor observed for 
some events. 



1 Introduction 



There are many evidences that faults and earthquakes interact, as suggested by calculations of stress 
redistribution [1], elastodynamic propagation of ruptures using laboratory-based friction law [2, 4], 
simplified models of multiple faults [5, 6], as well as general constraints of kinematic and geometric 
compatibility of the deformations [7]. Maybe the simplest mechanism for earthquake interaction 
involves stress re-distribution, both static [8, 1] and dynamical [9] associated with a given earth- 
quake modeled as a set of dislocations or cracks. In this simple mechanical view, earthquakes cast 
stress shadows in lobes of stress unloading [10, 8] and increase the probability of rupture in zones 
of stress increase [11], according to the laws of linear elasticity. These elastic stress transfer models 
are useful for their conceptual simplicity and are increasingly used. Notwithstanding their extended 
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use, the calculations of stress transfer have large uncertainties stemming from (i) the usually poorly 
known geometry of the rupture surfaces, (ii) the unconstrained homogeneity and amplitude of the 
stress drop and/or of the slip distribution on the fault plane, (iii) the use of simplified models of 
the crust (3D semi-infinite, or thin elastic plate, or plate coupled to a semi-infinite visco-elastic 
asthenosphere, etc.), and (iv) the unknown direction and amplitude of the absolute stress field that 
pre-existed before the event, including its possible spatial inhomogeneity. 

Such elastic stress transfer models seem unable to account for a growing phenomenology of 
long-range earthquake interactions. For instance, many large earthquakes have been preceded by 
an increase in the number of intermediate sized events over very broad areas [12, 13]. The relation 
between these intermediate sized events and the subsequent main event has only recently been rec- 
ognized on a large scale because the precursory events occur over such a large area that they do not 
fit prior definitions of foreshocks [14]. In particular, the 11 earthquakes in California with magni- 
tudes greater than 6.8 in the last century are associated with an increase of precursory intermediate 
magnitude earthquakes measured in a running time window of five years [15]. What is strange 
about the result is that the precursory pattern occured with distances of the order of 300 to 500 km 
from the futur epicenter, i.e. at distances up to ten times larger that the size of the futur earthquake 
rupture. Furthermore, the increased intermediate magnitude activity switched off rapidly after a 
big earthquake in about half of the cases. This implies that stress changes due to an earthquake of 
rupture dimension as small as 35 km can influence the stress distribution to distances more than 
ten times its size. These observations of earthquake-earthquake interactions over long times and 
large spatial separations have been strengthened by several other works on different catalogs using 
a variety of techniques [13, 16]. These results defy usual mechanical models of linear elasticity and 
one proposed explanation is that seismic cycles represent the approach to and retreat from a critical 
state of a fault network [16, 17]. Within the critical earthquake concept, the anomalous long-range 
interactions between earthquakes reflect the increasing stress-stress correlation length upon the ap- 
proach of the critical earthquake. Another explanation involves dynamical stress triggering [18] 
(see however [19]). Additional seismic, geophysical, and hydrogeological observations [20] cannot 
be accounted for by using models derived from the elastic stress transfer mechanism. In particu- 
lar, standard poro-elastic models underestimate grossly the observed amplitudes of hydrogeological 
anomalous rises and drops in wells at large distances from earthquakes. 

2 Mechanical Model of the Earth's Crust 

Here, we investigate the hypothesis, and its imphcations for the above observations, that the crust 
is a nonlinear elastic medium characterized by an asymmetric response to compressive versus ex- 
tensive perturbations around the lithostatic stress. We call this a "threshold nonlinearity." This 
nonlinearity stems from a mechanically-justified argumentation based on the fact that the Earth's 
crust at seismogenic depth is crisscrossed by joints, cracks or faults at many different scales filled 
with drained fluid in contact with delocahzed reservoirs at pressures close to the lithostatic pressure. 
It has been argued that rock permeability and thus microcracking adjusts itself, so that fluid pres- 
sure is always close to rock pressure irrespective of the extend of hydration/dehydration [21, 23]. 
One possible mechanism for this involves a time-dependent process that relates fluid pressure, flow 
pathways and fluid volumes [24]. 
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2.1 Presence and role of fluids 



Indeed, a lot of data collectively support the existence of significant fluid circulation to crustal 
depths of at least 10 — 15km. Much attention has been devoted to the role of overpressurized fluid 
[25, 26, 27, 28, 29, 30]. It is more and more recognized that fluids play an essential role in virtu- 
ally all crustal processes. Ref.[31] reviews the historical development of the conciousness among 
researchers of the ubiquitous presence and importance of fluids within the crust. Numerous exam- 
ples exist that demonstrate water as an active agent of the mechanical, chemical [32] and thermal 
processes that control many geologic processes that operate within the crust [21, 22]. The bulk of 
available information on the behavior of fluids comes from observations of exposed rocks that once 
resided at deeper crustal levels. In any case, present day surface exposed metamorphic rocks indi- 
cate that, at all crustal levels, fluids have been present in significant volume. Because the porosity 
of metamorphic rock is probably less than 1%, the high volume of calculated fluid necessary to pro- 
duce observed chemical changes suggests that fluid must have been replenished thousands of times. 
It has also been proposed that gold-quartz vein fields in metamorphic terranes provide evidence 
for the involvement of large volumes of fluids during faulting and may be the product of seismic 
processes [33]. Water is also released from transformed minerals. For instance, Montmorillonite 
changes to illite with a release of free water from the clay structure at approximately the same depth 
as the first occurrence of the anomalous pore pressure [34]. This is the most commonly discussed 
example of hydration and dehydration of minerals changing the fluid mass and the pore pressure. 
Fluids have been directly sampled at about 11 km by the Soviets at the Kola Peninsula drillhole. 

Many observations suggest that there are massive crustal fluid displacements correlated with 
seismic events. Among them, one can cite the fault-valve mechanism [35] or the migration and 
diffusion of aftershocks [36]. Several mechanisms have been proposed for the mechanical effect 
of fluids to decrease compressive lithostatic stresses: mantle-derived source of fluids can maintain 
overpressure within a leaky fault [27]; laboratory sliding experiments on granite show that the slid- 
ing resistance of shear planes can be significantly decreased by pore sealing and compaction which 
prevent the communication of fluids between the porous deforming shear zone and the surrounding 
material [28]. Physico-chemical processes such as mineral dehydration during metamorphism may 
provide large fluid abundance over large areas [37, 22]. Such fluid presence or migration implies 
that cracks may open and close even at seismogenic depth, justifying the relevance of asymetric 
nonlinear elasticity for such depths as we elaborate below. 

2.2 Physical mechanism for the threshold nonlinearity 

From a mechanical point of view, the threshold nonlinearity we invoke is different from the ubiq- 
uitous nonlinearity of rocks. In the later, nonlinearity becomes important only under large defor- 
mations. Such nonlinear elasticity of rocks is well-documented from its nonlinear wave signatures 
[39]. In contrast, the nonlinearity we invoke is revealed for almost arbitrary small perturbations by 
the difference between compressive versus extensive perturbations around a mean lithostatic stress 
field. Many crustal rocks have a Young's modulus depending on confining pressure, in particular 
with a Young's modulus in tension smaller than the Young's modulus in compression in a ratio from 
1/2 to 1/10 [38]. 

The highly damaged upper crust, as described above, would behave as a standard continuous 
elastic medium if almost all cracks are kept closed and are thus completely transparent to the ap- 
plied tectonic stress. In a totally dry crust, this would occur at depths larger than about Ikm. 
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Near the surface, cracks can open under sufficiently tensile stresses, so that the rock (if we neglect 
crack growth phenomena) will display low elastic modulii (what we will call later a soft state). 
Under compressive (or weakly tensile) stresses, those same cracks will close so that the rock will 
display elastic modulii close to those of uncracked material (the hard state). This simple compres- 
sive/tensile asymmetry changes strongly the behaviour of the rock: the rheology is still elastic (in 
the sense that it is reversible) but it is nonlinear (a change of sign of applied stress does not only 
changes the sign of strain, it also changes its modulus). In addition, if we now take account of the 
possible slow growth of cracks under tensile stress, the nature of the nonUnearity will not change, 
but its amplitude will. 

Let us now consider the influence of fluids whose presence are pervasive in the crust, as sum- 
marized in section 2.1. Let us assume that the crust is saturated with fluids and that the network of 
cracks is sufficiently dense so that it behaves as a drained medium. For the sake of simplicity, we 
will also assume that the crust is characterized by a homogeneous spatial distribution of cracks, so 
that the permeability is uniform and isotropic. We will also suppose that there exists a horizontal 
interface at depth Zgeai where permeabihty is 0. The last ingredient is that the part of the crust below 
Zscai is connected to a reservoir of fluid. Many indirect observations suggest indeed the presence of 
large sources of fluids [40, 41]. Then, above depth Zgeai, water will be at hydrostatic pressure Ph 
(lower than lithostatic pressure Pi) so that cracks will be able to close if the tensile applied tectonic 
stress is smaller than Pi — P^ in modulus (cracks remain closed if tectonic stress is compressive). 
Below Zseah the scenario is different. Water trapped in cracks is now at lithostatic pressure. If there 
is no applied tectonic stress, fluid pressure inside cracks compensates exactly the lithostatic pres- 
sure and the medium is exactly at the hard/soft boundary. The net stress acting on any crack's lips is 
0, so that the crack doesn't grow. If the applied tectonic stress is compressive, cracks will close and 
fluids are expelled towards the reservoir: the crust is in the hard phase. If the applied tectonic stress 
is tensile, cracks will open and the crust is in the soft phase. If a crack becomes unstable, pressure 
drops within it, so that it tends to close to re-establish initial fluid pressure. If a crack grows very 
slowly, pressure within the crack will stay constant with time. We will from now on neglect the 
possible slow crack growths and assume that fluid pressure is constant in each crack and that the 
crack network geometry does not show any evolution with time or applied stress. Thus, below Zgeah 
the applied stress threshold to get from the hard to the soft state is 0. This explains our terminology 
of a "threshold nonlinearity." 

The remaining of the paper will investigate stress redistributions associated with earthquakes 
occurring below Zgeai- The existence of Zgeai is difficult to prove for the real crust. It would cor- 
respond, for instance, to a depth where the crack network geometry changes abruptly, for instance 
at the boundary between the sedimentary basin and cristalline rocks. We can also propose an al- 
ternate model in which the vertical permeability below Zgcai is much lower than the horizontal one, 
and would obtain the same kind of soft/hard transition. Note that much more complicated scenarii 
involving the chemistry of both fluids and rock matrix could be taken into account [22], but we 
choosed to neglect them in this purely mechanical paper. 

We then propose a simple central-force spring model in two dimensions to study the behaviour 
of such an asymmetric nonlinear medium when subjected to an infinitesimal internal source of stress 
or strain. We show that the spatial structure of the stress transfer associated with such a perturbation 
evolves with the strength of the nonlinearity (defined as the ratio of the springs' stiffness in tensile 
and in compressive states). Those changes are quantified in terms of the symmetry of the resulting 
stress field and of the decay rate of the amplitude of the stress perturbation with distance from the 
source. 
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3 Numerical Model 



Solving theoretical problems of nonlinear elasticity proves to be very tough, even for the simple 
asymmetry of our threshold nonlinearity. We thus choosed to solve a couple of simple problems 
related to seismology using numerical modelling. Stress, strain and material rigidity being 2"'^ or 
Ath rank tensors, and as there is an obvious and complex feedback between strain and rigidity in our 
model, we choosed to use a simple spring model to illustrate the concept and major consequences 
of the threshold nonlinear elasticity. 

A plate of size L by L is discretized onto a regular grid of mesh size a. In the following, we 
choose L = 2000A;m and a = Wkm. Each elementary cell is defined by 4 nodes, each node 
being shared between 4 different neighbouring cells (except on the boundary of the plate). Figure 

1 shows the mechanical structure defined for each cell: each node is connected to its two nearest 
neighbours by springs of stiffness Ki (those springs indeed define the 4 edges of the cell). As those 
are central force springs, the shear modulus of such a cell is indeed 0. To get shear elasticity [42], 

2 springs are added along the cell's diagonals, such that each node is also connected to its next- 
nearest neighbours. The stiffness of those diagonal springs is K2. Once the plate is discretized with 
such cells, it can be shown that the plate behaves as an isotropic elastic medium if and only if we 
have K2 = Ki/2 [44]. The two independent elastic modulii of the plate can then be shown to be 
X = fi = Ki/2 for the two Lame coefficients, thus yielding E = (4/3)i^i = (8/3)^ for the Young 
modulus and = 1/3 for the Poisson coefficient (note that we are dealing with a pure 2D model). 
Of course, many other geometries of the elementary cell are possible, but we had to choose square 
cells which help to handle more easily with boundary conditions used in the problems we want to 
study. 

The relationships we just defined are true only if the springs are symmetric, i.e., their stiffnesses 
is the same under tensile or compressive states. The next and last step to model our nonlinear 
threshold elastic rheology is to impose that the stiffness of each spring can vary with its length. 
Thus, if a spring is shortened, its stiffness will be, say, K. If a spring is lengthened, its stiffness 
is lower and taken equal to aK, with a < 1. Each spring represents an effective volume filled 
with a uniform and isotropic distribution of cracks with sizes smaller than the representative mesh 
size. Real damage in the crust is of course much more complex, with anisotropic, space and scale 
dependence. These complications are neglected in our first exploration. We can obtain an order 
of magnitude estimate of the density of fluid-filled faults associated with a given asymmetric co- 
efficient a, using the effective medium calculations in [43]. To simplify, let us assume X = fi. 
Then, a = A^/Ai = 1/(1-1- 5d/2), where is the Lame coefficient of the damaged material and 
d = N{1/ Lf" is the density of faults (assumed identical) of radius ^ in a cube of volume L^. N 
is the number of faults in that volume. For instance, we need about 11 faults of size L/3 to get 
a = 0.5. Such estimate must however be taken with caution since the effective medium calculation 
a = 1/(1 + 5d/2) is valid only for small crack densities, while any piece of rock and the real 
crust are crisscrossed by many faults at many length scales, most of them being healed at varying 
degrees. We think that values of a significantly smaller than 1 should thus not be excluded. It is 
also probably that a is not uniform within the crust and can be expected to reflect the past history 
of deformations and ruptures. 

The ratio a of the extensive over compressive elastic coefficient can also be seen as equivalent 
to 1 — Djn in damage mechanics, where Dm is the scalar damage variable. If the spring isn't 
damaged, then D = 0, so that a = I and the stiffness is the same under tensile and compressive 
states. If the spring is totally damaged (near failure), D^ is close to 1, so that a ~ and the 
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stiffness of the spring in tensile state vanishes, while it is still K if the spring is compressed. Under 
arbitrary loading conditions, some springs in the plate will be in tensile state, while others will be 
in compressive state, so that it is difficult to analytically compute the stiffness tensor of the whole 
plate when a < 1. In addition, the stiffness tensor may feature more than 2 independent moduUi. 
This justifies the use of an iterative numerical method as described in the next section. 

4 Numerical Method 

The method we use belongs to the so-called iterative 'type-writer' methods. The first step consists in 
defining boundary conditions, i.e. to fix displacements and/or applied forces on set of nodes. As we 
are dealing with statics, such forces and displacements are held constants throughout the numerical 
process. We then consider, say, the node on the top left comer of the mesh and, according to 
the forces/displacements applied to this node and his nearest and next-nearest neighbours, we can 
compute the net force acting on that node. As we are dealing with a statics problem, the net force 
acting on that node at equilibrium should vanish. If the force vector we computed for that node is 
not , then we move it in the direction of the force vector to decrease the net force. We store the 
new position of this node and get to its right neighbour and follow the same scheme. We thus sweep 
the mesh line by line down to the node on the bottom right, and iterate the same operations again 
from the node on top left. As iterations accumulate, the net force acting on each node decreases, 
and we stop the process once all nodes are subjected to a force whose modulus is under a given 
threshold. This threshold is chosen such that the modulus of the incremental displacement necessary 
to decrease the net force modulus on each node is of the order of the accuracy of the computer, 
namely about 10~^^. To avoid any spurious result due to the direction of type- writing, iterations 
alternatively begin on each of the four corners of the plate on either horizontal or vertical directions, 
right or left, up or down. Once the 'type-writer' is stopped, displacements of all nodes relative to 
their initial positions are stored. Those positions also allow to compute, for a given node, the forces 
exerted on it by each of its neighbours. Those last quantities allow to compute the full stress tensor 
at that node [44], which is also stored. Computations of forces transmitted by a spring from one 
node to another take account of the spring stiffness which, as described in the previous section, 
depends on the state (stretched or shortened) of the spring through the value of a which is kept 
constant throughout the network. 

Starting from given boundary conditions, we solve a statics problem, which means that we 
do not take account of neither wave propagation nor fluid migration. We just compute the final 
equilibrium solution. This point will be discussed at the end of the paper when considering the 
application of this model to real Earth data. 

5 Earthquake Modeling 

In this paper, we want to examine the main differences between the stress field pattern generated by 
an earthquake in our nonlinear threshold elastic medium and in a standard linear elastic medium. 
The first step is thus to define mesh parameters representative of the real crust, the second one is to 
define what is an earthquake source in such a model. 
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5.1 Parametering the Earth's crust 

Our numerical model is strictly a 2D one, as dealing with the third dimension would lead to serious 
memory and computing time problems, and we would have to handle boundary problems such as 
the free surface and couphng with lower viscous layers. We can however choose mesh parame- 
ters such that, using relationships linking plane elasticity to "iD elasticity, we can model realistic 
crust properties (however neglecting boundary problems). The size of the plate as well of the cells 
has been previously given. We fixed the (virtual) thickness of the plate to be lO/cm, so that it 
corresponds roughly to the thickness of the seismogenic zone within the crust. We then choosed 
Kc = 5 ■ 10^^ Nm~^ (where Kc is the value of Ki when springs are compressed), so that the 3D 
elastic moduli! become E = 6.25 • 10^° Pa, A = /x = 2.5 • 10^° Pa, and u = 0.25, which are close 
to usual moduUi measured in rock mechanics experiments. The value of the assymmetry parameter 
a will be varied through several numerical experiments from 1 (standard isotropic elasticity) down 
to 0.01 (strong asymmetry of the elastic response under extension versus compression). 

5.2 Earthquake Source modeUing 

Earthquake source theory has until now been theoretically studied within the framework of linear 
elasticity. This allows one to use very powerful tools such as the representation theorem and Green 
functions. An earthquake can then be viewed equivalently as a displacement discontinuity across 
a fault plane, or a distribution of double-couples and dipoles of forces along the same plane in a 
continuous medium [3, 45, 46]. 

In the case of a fault of finite dimensions, if we assume that the stress drop is uniform along the 
fault, then we deal with a crack problem. If we assume that the displacement discontinuity across 
the fault is uniform, then we deal with a dislocation problem. At distances from the fault much 
larger than its size, and in the case of linear elasticity, both models yield the same spatial patterns 
of stress and displacement fields, which are linked by Hooke's law. This will be illustrated below 
in our diagrams obtained for the symmetric elasticity case a = 0. 

In the nonlinear case, it is easy to show that the representation theorem fails to apply, and then so 
does the Green function concept. This stems from the fact that the principle of linear superposition 
fails in the presence of nonlinearity. It follows that even the most simple earthquake source problem 
has to be defined either as a crack or a dislocation problem, and both problems should give different 
stress and displacement patterns at long wavelengths. We will, in this prehminary work, study only 
pointwise sources, which will allow simple comparisons with elementary solutions obtained from 
linear elasticity. 

We will thus study two cases: 

(i) an initially stress-free medium within which a single cell (located at its center) is subjected to 
the following stress field tensor: axx = (^yy = while axy 7^ (and is hereafter called point- 
wise shear stress load or crack model) - this model is reminiscent of the standard dynamical 
model of an event in standard hnear elasticity [3], 

(ii) an initially stress-free medium within which a single cell is subjected to a pure shear strain 
field: exx = f^yy = while exy 7^ (hereafter called pointwise shear strain load or disloca- 
tion problem) - this model rather views the event as a shear displacement discontinuity. 

In both cases, the corresponding infinitesimal planar defect (the source of the earthquake) suffers 
from undeterminacy and is oriented either along the Ox direction (plane Px) or along the Oy direc- 
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tion (plane Py). In the first case, the sUp discontinuity is dextral, and it is sinistral in the second 
case. 

5.3 Quantitative source parameters 

The small scale of our mechanical model is that of a cell, and this is thus the smallest scale we have 
to deal with in order to model an earthquake source. Figure 2 shows the source cell and 4 different 
vectors originating from each of its 4 comers. 

In the case of the pure shear stress load model, each vector represents a force applied to the 
corresponding node. All forces have the same modulus so that the stress tensor within the central 
cell indeed corresponds to the one we defined above. The same set of forces is applied whatever the 
value of a. 

In the case of the pure shear strain load model, each vector represents a displacement applied to 
the corresponding node. All displacements have the same modulus so that the strain tensor within 
the central cell indeed corresponds to the one we defined above. The same set of displacements is 
applied whatever the value of a. 

In order to ensure that we can compare results obtained from both types of boundary conditions, 
we have to fuUfil a very simple condition: the stress and displacement fields must be identical in the 
classic linear case, i.e. when a = 1. In the pure shear stress load model, we imposed the modulus 
of each applied force equal to F = 7.1 • 10^^ iV, so that is corresponds to an event of scalar moment 
Mo = 7.1 • 10^^ Nm, i.e., of magnitude ^ 6.5. We computed the displacement field for a = 1, and 
observed that the magnitude of the displacement at each node of the source cell was 0.526640588m. 
To be perfectly consistent, we thus impose this displacement amplitude at each source cell node in 
the case of the pure shear strain load model. 

6 Displacement field at the source 

In the pointwise dislocation case, displacements are held constant at the source whatever the value 
of a. In the pointwise crack model, only forces are kept constant as a is changed, and displacements 
are expected to vary as the asymmetry increases (i.e. a decreases). We have already pointed out 
that both boundary conditions assume that the mechanical defect is either parallel to Px or to Py. 
Indeed, we will show that the displacement (hereafter named ?/„) normal to each of these conjugate 
defects is of the same type (opening) while the shear displacement (hereafter Us) along each of them 
is of different type: dextral along P^ and sinistral along Py. That said, we will focus only on the 
modulii of those displacements. 

How do we obtain m„ and Us? Indeed, according to the orientations of the conjugate plane 
defects, we have to compute displacements at points A, B, C or D, whereas our model provides 
solutions at nodes 1,2,3 and 4 (see Figure 2). Displacements at points Aio D thus com- 
puted through bilinear interpolations within the cell. We then define u„ = Ux{B) = —Ux{D) = 
—Uy{C) = Uy{A) and Ug = Ux{A) = Uy{B) = —Ux{C) = —Uy{D) as a result of the symmetries 
of the system. 

Figure 3 shows the variations of both Un and Us with a (in fact, it shows the values of displace- 
ments discontinuities across the crack, i.e. 2un and 2us). For a = 1, we find that Un is very close 
to 0, so that the displacement along the pointwise crack lips is pure shear. As a decreases from 1 to 
close to 0, the shear displacements increase by a factor of about 4.5. This is perfectly understand- 
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able, as some springs are soUicited in tension, leading to a decrease in their stiffness. This decrease, 
in the presence of constant forces, implies that displacements increase. 

More surprising is the behaviour of normal displacements, which increase drastically as a de- 
creases, tending to be about 2.7m when a tends to (i.e. about half the shear displacements). 
For a = 0.1 we have Un/ug close to 1/3. Moreover, and Ug are both positive, which signifies 
that, under the shear stress load assumption, the defect opens when the medium is asymmetric. In 
seismological words, this means that the static moment tensor of the source has a non-vanishing 
trace and can thus be decomposed into an isotropic part and a deviatoric one. Despite the ob- 
servation that most earthquake sources are thought to be well modelled by the deviatoric part 
alone, a few catalogues report isotropic components. Several mechanisms have been invoked to 
explain a non-vanishing isotropic component of the seismic moment. The standard explanation for 
non-double-couple components relies on the fault zone irregularity [47]. Some earthquakes with 
non-double-couple mechanisms have been claimed not to be explained solely by such a composite 
rupture [48, 22]. It is then important to note that the seismic moment tensor reported in catalogues 
is a variable that quantifies static displacements at the source. A static dilational strain at the source 
can occur even when the dynamic representation of the source is a pure double couple of forces 
(yielding a stress tensor with zero trace). 

Figure 3 also reports the shear strain that can be measured at the source cell as well as the 
relative dilation of its surface dS/ S. Both quantities of course corroborate the previous results on 
Us and Un- 

7 The stress field 

We will now focus on the structure of the stress field generated by our pointwise earthquakes. The 
stress field will be studied at large wavelengths, i.e. at distances from the source cell of more than a 
few cell sizes. This gives the rate of stress decay with distance from the source and thus the range 
of interactions between events. We will in the following implicitly assume that the plate is affected 
by many other faults that are locked and oriented in the same direction (say, along direction with 
potential dextral displacement along the fault). The source cell is one of such fault producing an 
event. We then study the effect of this event on all other faults in the plate. 

7.1 Spatial patterns 

Figures 4 to 6 show the variation of the shear stress component axy within the plate near the source 
cell. A positive variation signifies that this stress component increases. On each figure, the panels on 

the left represent patterns obtained with the pure shear stress load hypothesis, while the panels on the 
right represent patterns obtained with the pure shear strain load hypothesis. Each row corresponds 
to a different value of a, i.e., to a different degree of elastic asymmetry between extensive and 
compressive deformation. 

In the case a = 1, both boundary conditions yield exactly the same spatial pattern, as expected. 
This is in agreement with the fact that, in linear elasticity, both boundary conditions are equivalent. 
We obtain 8 lobes of identical shapes within which stress amplitude alternates from positive to 
negative (red and blue lobes, repectively). 

As a decreases, the symmetry of the patterns decreases: some lobes are rotating, some are 
widening while others are narrowing. We will study stress variation in those lobes in a subsequent 
section. The most striking observation is however that the spatial structure of the patterns is the 
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same for both types of boundary conditions for a given value of a, notwithstanding a significantly 
smaller stress amplitude in the pure shear strain load boundary condition. 

Figures 7 to 9 show the variation of the stress component axx within the plate. The same kind 
of comments apply here as for Gxy This is also the case for component dyy which is shown in 
Figures 10 to 12. 

We have thus introduced a mechanical asymmetry at the microscopical level (i.e. at the spring 
scale), for all springs' orientations corresponding to an isotropic asymmetry, which translates into a 
loss of symmetry at the macroscopical scale. We shall quantify this loss of symmetry more precisely 
in the next section. 

7.2 Decay of the stress amplitude away from the source 

If we consider the center of the source cell as the origin of our frame, every point within the plate 
can be located using polar coordinates (r, 9), where r is the distance to the origin (the source cell 
center), and 9 is an azimuth measured clockwise from the Ox axis. We can then, for a fixed 6, look 
at the decay rate with r of the modulus of any of the stress tensor components. To avoid problems 
due to the finite size of the source and of the whole plate, we quantify the decay of the modulus of 
any stress component by a power-law of the type r~'^, within a distance interval bracketed within 
a few cell to a few tens of cell sizes. We then repeat the same computation for different values of 
9. We then change the value of the asymmetry factor a and obtain the corresponding dependences 
7(^) for each value of a. 

Figures 13 to 15 show the variation of 7 with 0, for different values of a, quantifying the decay 
of the Gxy component. Each frame features two curves, one corresponding to the pure shear stress 
condition at the source (in red), the other one to the pure shear strain condition (in blue). 

When a = 1, 7 is very close to 2, which is the theoretical value predicted by planar elasticity. 
There are some small fluctuations around that value, the largest ones being obtained for values of the 
azimuth corresponding to a change of sign of the stress, i.e. where the stress itself almost vanishes. 
In those peculiar directions, the determination of the exponent is very unstable. 

When a decreases, 7 values can reach values very different from 2. Some of those values 
correspond to azimuths where stress vanishes (and are thus spurious), but others reflect genuine 
consequences of the nonlinearity of the medium. One can see that exponents can thus reach values 
larger than 2 (reflecting a very rapid decay and thus short interaction range), but that they can also 
get down to values around or lower than 1, leading to very large interaction ranges. When a = 0.01, 
one should not consider negative values of 7 too seriously as such negative values would imply that 
the stress increases with distance. The increase is occurring only over a finite distance range and 
gives way to a decrease at larger distance. The measured exponent is thus only vahd at very short 
distances and is not an asymptotic value. We are unable for those cases to quantify accurately the 
value of the asymptotic 7 due to finite size effects. For that a values, one should thus consider that 
the asymptotic value of 7 varies en general between and 1. 

The found values of 7 as a function of 9 also reveal that exponents do not vary significantly 
with the conditions imposed at the source, which constitutes another surprise. However, stresses 
obtained in the pure shear strain condition are lower than in the pure shear stress condition, as the 
prefactor of the power-law decay is found smaller than for the pure shear stress case. 

Figures 16 to 18, as well as Figures 19 to 21 show the same results for components axx and 
Gyy. For both components, and for a = 1, we recover the theoretical value 7 = 2 for any 9. As 
a decreases, the exponents can take very different values, including some which imply very long 
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range decay. We observe again that the exponents do not vary with the type of loading at the source. 
We also checked that the explaination for negative values of exponents was the same as for Uxy 
Thus, for low a values, 7 decreases to values between and 1. 

8 Discussion 

This idea of mechanical asymmetry, and/or of the feedback between local damage and stress decay 
from perturbative sources is not new but, to the best of our knowledge, it is the first time that it is 
implemented in a real 2D plane elastic problem applied to Earth mechanics and seismotectonics. 
For example, [51] gives the analytical solution for the stress field and for the dependence 7(0;) in 
a nonUnear asymmetric elastic medium in the case of antiplane mode m loading, which is thus 
the scalar equivalent to the problem studied here. In the antiplane case, there is only one stress 
component and a single exponent 7(a). In this antiplane case, it can be shown analytically that the 
exponent 7(a) is indeed decreasing from the value 2 for a = 1 to smaller values as a decreases. 
But there is not dependence on azimuth for this scalar case. 

The existence of an asymmtry in the crust elasticity has been proposed on the basis of observa- 
tions of the Manyi {AIw = 7.6) earthquake [49]. Using SAR interferometry data, Peltzer et al. [49] 
interpreted the mismatch between the displacement across each side of the left-lateral strike-slip 
fault as due to a mechanical asymmetry between dilational and compressional quadrants. Using 
a first-order perturbative calculation, they estimated a coefficient 1/4 < a < 1/2 to explain the 
observed displacement asymmetry. However, they did not consider the possibility the asymmetry 
could modify the long range decay rate of the stress field. This in turn can modify the future seismic 
history in the neighborhood of this event. Their computation showed that the asymmetric effect was 
probably confined in the very shallow part of the crust which, if true, implies that the stress transfer 
at seismogenic depths after this event obeys standard linear elastic solutions. 

The fact that this model should be relevant for shallow crustal mechanics is rather obvious (as 
shown from the previous field example as well as from the short discussion at the beginning of 
this paper). The important question is to check if this model also holds (at least in limited spatial 
domains) at depth. In that case, stress transfer case-studies should take account of this asymmetry 
effect, which can greatly enhance the distance at which a given event can trigger another one. 
Testing this hypothesis is not simple, as the rheology we assumed is nonlinear, which means that 
the effect of successive events can not be simply added. Stress field evolutions with time may 
then have much sharper transitions in space and time than predicted by models involving linear 
elasticity, a behaviour reminiscent of the mechanics of granular media [50]. The consequence is 
that, to compare with the standard stress transfer mechanism [ 1 ] , we need to know in details the state 
of stress within the crust prior to an event to map predicted stress transfer lobes onto aftershocks 
location catalogs. 

Another possibihty for testing our hypothesis would be to study the statistics of seismic moment 
tensors of events, as we saw that, if we assume that the source can be describe dynamically as a 
pure double couple of forces, this tensor should display a non-vanishing isotropic part. However, 
the seismic moment tensor is computed from seismic wave observations (i.e. of dynamical nature), 
and not from static displacements in situ at depth. Morover, the time at which the solutions we 
computed really hold depends on the diffusion properties of fluids in rock. This is why such a 
way of testing would indeed imply to compute the whole time-dependent dynamical asymmetric 
poro-elastic solution to really propose quantitative results allowing a reliable comparison. 
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Other data that could be used to test the model are the SAR data, in the spirit of the work of 
[49]. Interpretation of SAR data done after a large event could prove the pertinence of our model, 
provided that stress field evolution after the event is not modified by other events or by the far-field 

loading of the plate. 

Proving or refuting this model is of prime importance for the understanding of the spatio- 
temporal patterns of earthquakes, including those preceeding a large and potentially destructive 
event. We already discussed the fact that it could improve the potential of prediction methods based 
on concepts such as stress transfer. But it has even deeper implications on the hope of predicting 
large events from the behavior of the statistics of populations of shocks preceeding that event. In 
another numerical work, Ref. [52] studied the progressive damage of a fault plane before its macro- 
scopic rupture. Their model employs cellular automaton techniques to simulate tectonic loading, 
rupture events and strain redistribution. Note that in that case, strain is equivalent to stress. The 
elastodynamic Green function for stress/strain redistribution is taken to vary as l/r^, where p is a 
parameter which is varied. The systems displays two different regimes depending on the p value. 
For p < 2, large events are preceeded by a clear power-law acceleration of energy release of the 
system, together with the growth of strain energy correlations. This is of course reminiscent of the 
critical earthquake hypothesis [13, 16, 17]. For p larger than 2, the trend of energy release before a 
large event is linear. This means that the lower p is, the more predictable is the large event, using 
time series of precursory energy release. Their model does not map exactly to ours,, but we could 
expect that if a is under a certain threshold (still to be determined), then 7 would be low enough 
for the critical earthquake scenario to apply, making large events predictable from time series. In 
the other hand, if in some areas a is above that threshold, then the local tectonic domain would 
belong to the other regime, and large events would be unpredictable. If the predictability of large 
events relies, as suggested by [52], on the exponent of the Green function of stress transfer, then 
we have pointed out a very simple physical mechanism allowing to tune that exponent. Large scale 
fluctuations of fluid pressure from lithostatic to infra-lithostatic could then explain why in some 
cases large events are preceeded by strain energy release acceleration, while the opposite holds in 
other cases. 
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Figure 1 : Mechanical structure of an elementary cell of the model. Each cell is defined by 4 nodes 
linked by springs along edges and both diagonals. Springs on the edges are of stiffness Ki, while 
diagonal springs are of stiffness K2 = Ki/2 (see text). 
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Figure 2: Representation of the elementary source in a unit cell. Vectors represent either forces 
of equal modulii (crack model) or displacement of same amplitude (dislocation model). Nodes of 
the cell are labelled 1, 2, 3, 4. Dotted lines represent the virtual plane defect along which force or 
displacement discontinuities are imposed to model an earthquake. Those discontinuities have to be 
computed at points labelled A, B, C, D (see text). 
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Figure 3: Left panel: shear (circles) and normal (triangles) displacements discontinuities along 
planar defects within the source cell as a function of a in semilogscale. Results for a = are not 
shown as the normal discontinuity vanishes. Right panel: volumetric (circles) and shear (triangles) 
strains of the source cell as a function of a in semilogscale (a is the ratio of the elastic modulus 
in extension over the elastic modulus in compression). Results for a = are not shown as the 
volumetric strain vanishes in this case. All results are obtained for the pure shear stress loading 
case. 
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Figure 4: Map of the shear stress transfer axy for a = 1 and a = 0.8 in the pure shear stress load 
case (left panels) and in the pure shear strain loading case (right panels). The source is located at 
(0, 0) and all space units are kilometers. 



18 




-100 -50 50 -100 -50 50 



Figure 5: Same as Fig. 4 for a = 0.6 and a = 0.4. 
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Figure 7: Map of the stress transfer Oxx for a = 1 and a = 0.8 in the pure shear stress load case 
(left panels) and in the pure shear strain loading case (right panels). The source is located in (0, 0) 
and all space units are kilometers. 
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Figure 8: Same as Fig. 7 for a = 0.6 and a = 0.4. 
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Figure 9: Same as Fig. 7 for a = 0.2 and a = 0.01. 
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Figure 10: Map of the stress transfer ayy for a = 1 and a = 0.8 in the pure shear stress load case 
(left panels) and in the pure shear strain loading case (right panels). The source is located in (0, 0) 
and all space units are kilometers. 
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Figure 11: Same as Fig. 10 for a = 0.6 and a = 0.4. 
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Figure 13: Variation of the decay exponent 7 with azimuth 9 when a = 1 measured for the stress 
component a-xy The values of 7 obtained with pure shear stress loading and pure shear strain 
loading superimpose exactly. All values of 7 are close to the theoretical value 7 = 2 for symmetric 
elasticity, except at stress lobes borders where it is ill-defined numerically and exhibits spurious 
peak values. 
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Figure 14: Variation of the decay exponent 7 with azimuth ^ for a = 0.5 for the stress component 
axy. The red (respectively blue) curve is for pure shear stress (respectively pure shear strain) load- 
ing. Peak values are spurious and correspond to lobes boundaries. The exponents are essentially 
identical for both types of loading conditions. 
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Figure 16: Variation of the decay exponent 7 with azimuth for a = 1 for the stress component a^x- 
The values of 7 obtained with pure shear stress and pure shear strain loading superimpose exactly. 
7 is found close to the theoretical value 2, except at stress lobe borders where it is ill-defined and 
exhibits spurious peak values. 
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Figure 17: Same as Fig. 16 for a = 0.5. The red (respectively blue) curve is for pure shear stress 
(respectively strain) loading. Peak values are spurious and correspond to lobe boundaries. The 
exponents are very similar for both types of loading conditions. 
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Figure 18: Same as Fig. 16 for a = 0.01. 
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Figure 19: Variation of the decay exponent 7 with azimuth 9 for a = 1 for the stress component 
ayy. The values of 7 obtained with pure shear stress and pure shear strain loading superimpose 
exactly. All values of 7 are close to the theoretical value 7 = 2, except at stress lobe borders where 
they are ill-defined and exhibit spurious peak values. 
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Figure 20: Same as Fig. 19 for a = 0.5. The red (respectively blue) curve is for pure shear stress 
(respectively strain) loading. Peak values are again spurious and correspond to lobe boundaries. 
The exponents are very similar for both types of loading conditions. 



34 



3 




100 



Figure 21: Same as Fig. 20 for a = 0.01. 
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